#### Meta Data ####

# article: "The Ambivalent Effect of Autocratization on Domestic Terrorism" 
# journal: "Studies in Conflict & Terrorism"
# authors: "Lott, Lars; Croissant, Aurel; Trinn. Christoph"
# date: 2023-10-02
# written under "R version 4.3.1 (2023-06-16 ucrt)"

#### Preliminaries ####
R.version$version.string

# clear workspace
rm(list=ls())

# set working directory
getwd()

#### Preliminaries ####

library(MASS)
library(multiwayvcov)
library(ggeffects)
library(tidyverse)
library(ggpubr)
library(here)
library(haven)
library(imputeTS)
library(modelsummary)
library(data.table)
library(patchwork)

#### functions ####

panel_match_path <- "PanelMatch/R/"

source(file.path(panel_match_path, "data.R"))
source(file.path(panel_match_path, "DisplayTreatment.R"))
source(file.path(panel_match_path, "exact_match.R"))
source(file.path(panel_match_path, "listwise_delete.R"))
source(file.path(panel_match_path, "matched_set_obj.R"))
source(file.path(panel_match_path, "matched_set_R.R"))
source(file.path(panel_match_path, "network_and_caliper.R"))
source(file.path(panel_match_path, "PanelEstimate.R"))
source(file.path(panel_match_path, "PanelEstimateObject.R"))
source(file.path(panel_match_path, "PanelMatch.R"))
source(file.path(panel_match_path, "PanelMatch-Package.R"))
source(file.path(panel_match_path, "PE_helpers.R"))
source(file.path(panel_match_path, "PE_lower_level.R"))
source(file.path(panel_match_path, "pm_helpers_r.R"))
source(file.path(panel_match_path, "RcppExports.R"))
source(file.path(panel_match_path, "utilities.R"))

#### Load Dataset ####

vdem <- readRDS("data/vdem12.rds") 

#### Define Treatment ####

vdem <- vdem %>%
  mutate(autocratization_treatment = aut_ep,
         democratic_breakdown = ifelse(aut_ep_outcome_agg == 1, 1, 0), 
         democratic_recession = ifelse(aut_ep_outcome_agg == 2, 1, 0), 
         regressed_autocracy = ifelse(aut_ep_outcome_agg == 3, 1, 0))

table(vdem$autocratization_treatment)
table(vdem$democratic_breakdown)
table(vdem$democratic_recession)
table(vdem$regressed_autocracy)

vdem <- vdem %>%
  mutate(democratic_breakdown_recession = ifelse(democratic_breakdown == 1 |democratic_recession == 1, 1,0 ))

table(vdem$democratic_breakdown_recession)



vdem <- vdem %>%
  mutate(log_count_attacks_dom = log(count_attacks_dom +0.01),
         log_count_attacks_int = log(count_attacks_int +0.01), 
         log_count_attacks = log(count_attacks +0.01),
         log_sum_nkill_dom = log(sum_nkill_dom +0.01),
         log_sum_nkill_int = log(sum_nkill_int +0.01),
         log_sum_nkill = log(sum_nkill +0.01)) %>%
  filter(year >= 1966)


vdem <- vdem %>%
  group_by(country_id) %>%
  mutate(log_count_attacks_dom = na_interpolation(log_count_attacks_dom, option = "spline",  maxgap = 1), 
         log_sum_nkill_dom = na_interpolation(log_sum_nkill_dom, option = "spline",  maxgap = 1), 
         count_attacks_dom = na_interpolation(count_attacks_dom, option = "spline",  maxgap = 1), 
         sum_nkill_dom = na_interpolation(sum_nkill_dom, option = "spline",  maxgap = 1), 
         log_count_attacks_dom = ifelse(log_count_attacks_dom<0, 0, log_count_attacks_dom),
         log_sum_nkill_dom = ifelse(log_sum_nkill_dom<0, 0, log_sum_nkill_dom),
         count_attacks_dom = ifelse(count_attacks_dom<0, 0, count_attacks_dom),
         sum_nkill_dom = ifelse(sum_nkill_dom<0, 0, sum_nkill_dom))

vdem <- vdem %>%
  dplyr::select(country_id, country_name, year, aut_ep, 
                aut_ep_id,  log_count_attacks_dom, log_count_attacks_int, log_count_attacks, 
                log_sum_nkill_dom, log_sum_nkill_int, log_sum_nkill, count_attacks_dom, count_attacks_int, count_attacks, 
                sum_nkill_dom, sum_nkill_int, sum_nkill,e_gdppc, e_pop,
                v2x_polyarchy, milper_extrapol, eprratio_extrapol, democratic_recession, democratic_breakdown, 
                civil_war_ucdp, v2x_regime, aut_ep_outcome, aut_ep_outcome_agg, 
                v2caviol, dem_ep, dem_ep_id, dem_ep_outcome_agg, dem_ep_outcome, e_regionpol, e_regionpol_6C, 
                autocratization_treatment, regressed_autocracy, democratic_breakdown_recession, ma_count_attacks_dom, ma_sum_nkill_dom, democracy_stock, 
                v2x_veracc, v2x_diagacc, v2x_horacc, e_pelifeex, e_peaveduc, e_area) %>%
  filter(year <=2020) %>%
  group_by(country_id) %>%
  fill(e_area)

vdem <- vdem %>%
  drop_na(autocratization_treatment) %>%
  mutate(e_gdppc_sqaured = e_gdppc^2, 
         e_area_log = log10(e_area))

summary(vdem$e_area_log)
summary(vdem$e_gdppc_sqaured)


#### Descriptive Statistics Table A1 ####

vdem_summary_data <- vdem %>%
  filter(year >=1970 & year <=2020) %>%
  dplyr::select(country_id, year, count_attacks_dom, 
                autocratization_treatment, democratic_breakdown_recession, regressed_autocracy,
                e_gdppc, milper_extrapol, eprratio_extrapol, 
                e_pop, democracy_stock, e_area_log, ma_count_attacks_dom)

datasummary_skim(vdem_summary_data, output = 'outputs/Table_A1.docx', fmt="%.4f")
datasummary_skim(vdem_summary_data, fmt="%.4f")

datasummary_balance(data = vdem_summary_data, ~ autocratization_treatment,  fmt="%.4f")
datasummary_balance(data = vdem_summary_data, ~ democratic_breakdown_recession, fmt="%.4f")
datasummary_balance(data = vdem_summary_data, ~ regressed_autocracy, fmt="%.4f")


#------------------------------------------------------------------------------#
#---------------------------------- Figure 1 ----------------------------------#
#------------------------------------------------------------------------------#

vdem_year_data <- vdem %>%
  group_by(year) %>%
  summarize(sum_aut_ep = (sum(aut_ep, na.rm = T)), 
            sum_terrorist_attacks = sum(count_attacks_dom, na.rm = T)) %>%
  filter(year >= 1970)

f1 <- ggplot(data = vdem_year_data, aes(x = year, y = sum_aut_ep)) +
  geom_line() +
  tidyquant::geom_ma(ma_fun = SMA, n = 5, color = "red") +      
  scale_y_continuous(breaks = seq(0, 10, by = 1)) +
  labs(x = "Year", y = "Number of Ongoing\nAutocratization Waves") +
  theme(legend.position = "none", 
        plot.title = element_text(hjust = 0.5)) +
  theme_pubr()

f2 <- ggplot(data = vdem_year_data, aes(x = year, y = sum_terrorist_attacks)) +
  geom_line(linetype = "twodash") +
  labs(x = "Year", y = "Number of Domestic\nTerrorist Attacks") +
  theme(legend.position = "none", 
        plot.title = element_text(hjust = 0.5)) +
  theme_pubr()

F1_main <- f1 / f2

ggsave(plot = F1_main, "outputs/Figure_1.png", dpi = 600, units = "cm", height = 12, width = 19 )


#------------------------------------------------------------------------------#
#---------------------------------- Figure 2 ----------------------------------#
#------------------------------------------------------------------------------#

# not plotted in R. We used Latex for this Figure


#------------------------------------------------------------------------------#
#---------------------------------- Figure 3 ----------------------------------#
#------------------------------------------------------------------------------#
vdem <- distinct(vdem, country_id, year, .keep_all= TRUE) # Distinct observations

vdem$year <- as.integer(vdem$year)


vdem <- as.data.frame(vdem)
class(vdem)

# Define Lags and leads for all

lag_2 <- c(1:2)
lead_5 <- c(0:5)


vdem$country_id <- as.integer(vdem$country_id)

summary(vdem$democratic_breakdown_recession)
table(vdem$democratic_breakdown_recession)


#### 1. Effect of democratic_breakdown_delta on Terrorism in all regimes #####

## Before Refinement ##
PM.results.none <- PanelMatch(lag = 2, time.id = "year", unit.id = "country_id", 
                              treatment = "democratic_breakdown_recession", refinement.method = "none", 
                              data = vdem, match.missing = FALSE, 
                              covs.formula = ~ I(lag(ma_count_attacks_dom, lag_2)), 
                              size.match = 10, qoi = "att" ,outcome.var = "count_attacks_dom",
                              lead = lead_5, forbid.treatment.reversal = FALSE, 
                              use.diagonal.variance.matrix = TRUE)

## Mahanlanobis Distance Matching ##
PM.results_mahalanobis <- PanelMatch(lag = 2, time.id = "year", unit.id = "country_id", 
                                     treatment = "democratic_breakdown_recession", refinement.method = "mahalanobis", 
                                     data = vdem, match.missing = FALSE, 
                                     covs.formula = ~ I(lag(e_gdppc, lag_2)) + 
                                       I(lag(milper_extrapol, lag_2)) + I(lag(eprratio_extrapol, lag_2)) + 
                                       I(lag(e_pop, lag_2)) +  I(lag(democracy_stock, lag_2)) + I(lag(e_area_log, lag_2)) +  
                                       I(lag(ma_count_attacks_dom, lag_2)), 
                                     size.match = 10, qoi = "att" ,outcome.var = "count_attacks_dom",
                                     lead = lead_5, forbid.treatment.reversal = FALSE, 
                                     use.diagonal.variance.matrix = TRUE)

## Propensity Score Matching ##
PM.results_psmatch <- PanelMatch(lag = 2, time.id = "year", unit.id = "country_id", 
                                 treatment = "democratic_breakdown_recession", refinement.method = "ps.match", 
                                 data = vdem, match.missing = FALSE, 
                                 covs.formula = ~ I(lag(e_gdppc, lag_2)) + 
                                   I(lag(milper_extrapol, lag_2)) + I(lag(eprratio_extrapol, lag_2)) + 
                                   I(lag(e_pop, lag_2)) +  I(lag(democracy_stock, lag_2)) + I(lag(e_area_log, lag_2)) +  
                                   I(lag(ma_count_attacks_dom, lag_2)), 
                                 size.match = 10, qoi = "att" ,outcome.var = "count_attacks_dom",
                                 lead = lead_5, forbid.treatment.reversal = FALSE,
                                 use.diagonal.variance.matrix = TRUE)

## Propensity Score Weighting ##
PM.results_psweight <- PanelMatch(lag = 2, time.id = "year", unit.id = "country_id", 
                                  treatment = "democratic_breakdown_recession", refinement.method = "ps.weight", 
                                  data = vdem, match.missing = FALSE, 
                                  covs.formula = ~ I(lag(e_gdppc, lag_2)) + 
                                    I(lag(milper_extrapol, lag_2)) + I(lag(eprratio_extrapol, lag_2)) + 
                                    I(lag(e_pop, lag_2)) +  I(lag(democracy_stock, lag_2)) + I(lag(e_area_log, lag_2)) +  
                                    I(lag(ma_count_attacks_dom, lag_2)), 
                                  size.match = 10, qoi = "att" ,outcome.var = "count_attacks_dom",
                                  lead = lead_5, forbid.treatment.reversal = FALSE,
                                  use.diagonal.variance.matrix = TRUE)

class(PM.results_psweight)

msets <- PM.results_psweight$att

print(msets) # 32 democratic breakdown

######################################################

## Display different matched sets ##
mset <- PM.results_psweight$att[10]

DisplayTreatment(unit.id = "country_id",
                 time.id = "year", legend.position = "none",
                 xlab = "year", ylab = "Country Code",
                 treatment = "democratic_breakdown_recession", data = vdem,
                 matched.set = mset, # this way we highlight the particular set
                 show.set.only = TRUE)

summary(PM.results_mahalanobis$att)
summary(PM.results_psweight$att)


par(mfrow=c(1,1)) 
plot(PM.results_psweight$att, main ="Distribution of Matched Set Sizes")

dev.copy(pdf,'results/Figure_3/Figure_Frequency_Distribution.pdf', width = 10, height = 5)
dev.off()


#### Check Improved Covariate Balance due to Matching over the Pre-Treatment Time ####

## Get covariate balance ##

df_none <- get_covariate_balance(PM.results.none$att,
                                 data = vdem,
                                 covariates = c("e_gdppc",  "milper_extrapol", "eprratio_extrapol",
                                                "e_pop", "ma_count_attacks_dom", "democracy_stock", "e_area_log"),
                                 plot = FALSE)

df_mahalanobis <- get_covariate_balance(PM.results_mahalanobis$att,
                                        data = vdem,
                                        covariates = c("e_gdppc",  "milper_extrapol", "eprratio_extrapol",
                                                       "e_pop", "ma_count_attacks_dom", "democracy_stock", "e_area_log"),
                                        plot = FALSE)


df_psmatch <- get_covariate_balance(PM.results_psmatch$att,
                                    data = vdem,
                                    covariates = c("e_gdppc",  "milper_extrapol", "eprratio_extrapol",
                                                   "e_pop", "ma_count_attacks_dom", "democracy_stock", "e_area_log"),
                                    plot = FALSE)


df_psweight <- get_covariate_balance(PM.results_psweight$att,
                                     data = vdem,
                                     covariates = c("e_gdppc",  "milper_extrapol", "eprratio_extrapol",
                                                    "e_pop", "ma_count_attacks_dom", "democracy_stock", "e_area_log"),
                                     plot = FALSE)


#Tidying outputs #
df_none <- as.data.frame(df_none)

df_none <- df_none %>%
  rownames_to_column("year") %>%
  pivot_longer(-year, names_to = "covariate", values_to = "sd_covariate") %>%
  mutate(year = case_when(year =="t_4" ~ -4, 
                          year =="t_3" ~ -3,
                          year =="t_2" ~ -2,
                          year =="t_1" ~ -1,
                          year == "t_0" ~ 0))

F1 <- ggplot(df_none, aes(x = year, y = sd_covariate, color = covariate)) +
  geom_line() +
  geom_hline(yintercept=0, color = "black") +
  labs(x = "", y = "", 
       title = "Before Refinement") +
  ylim(-0.8, 0.8) +
  theme(legend.position = "none", 
        plot.title = element_text(hjust = 0.5)) +
  theme_pubr()

df_mahalanobis <- as.data.frame(df_mahalanobis)

df_mahalanobis <- df_mahalanobis %>%
  rownames_to_column("year") %>%
  pivot_longer(-year, names_to = "covariate", values_to = "sd_covariate") %>%
  mutate(year = case_when(year =="t_4" ~ -4, 
                          year =="t_3" ~ -3,
                          year =="t_2" ~ -2,
                          year =="t_1" ~ -1,
                          year == "t_0" ~ 0))

F2 <- ggplot(df_mahalanobis, aes(x = year, y = sd_covariate, color = covariate)) +
  geom_line() +
  geom_hline(yintercept=0, color = "black") +
  labs(x = "", y = "", 
       title = "Mahanlanobis Distance\n Matching") +
  ylim(-0.8, 0.8) +
  theme(legend.position = "none", 
        plot.title = element_text(hjust = 0.5)) +
  theme_pubr()

df_psmatch <- as.data.frame(df_psmatch)

df_psmatch <- df_psmatch %>%
  rownames_to_column("year") %>%
  pivot_longer(-year, names_to = "covariate", values_to = "sd_covariate") %>%
  mutate(year = case_when(year =="t_4" ~ -4, 
                          year =="t_3" ~ -3,
                          year =="t_2" ~ -2,
                          year =="t_1" ~ -1,
                          year == "t_0" ~ 0))

F3 <- ggplot(df_psmatch, aes(x = year, y = sd_covariate, color = covariate)) +
  geom_line() +
  geom_hline(yintercept=0, color = "black") +
  labs(x = "", y = "", 
       title = "Propensity Score\n Matching") +
  ylim(-0.8, 0.8) +
  theme(legend.position = "none", 
        plot.title = element_text(hjust = 0.5)) +
  theme_pubr()


df_psweight <- as.data.frame(df_psweight)

df_psweight <- df_psweight %>%
  rownames_to_column("year") %>%
  pivot_longer(-year, names_to = "covariate", values_to = "sd_covariate") %>%
  mutate(year = case_when(year =="t_4" ~ -4, 
                          year =="t_3" ~ -3,
                          year =="t_2" ~ -2,
                          year =="t_1" ~ -1,
                          year == "t_0" ~ 0))

F4 <- ggplot(df_psweight, aes(x = year, y = sd_covariate, color = covariate)) +
  geom_line() +
  geom_hline(yintercept=0, color = "black") +
  labs(x = "", y = "", 
       title = "Propensity Score\n Weighting") +
  ylim(-0.8, 0.8) +
  theme(legend.position = "none", 
        plot.title = element_text(hjust = 0.5)) +
  theme_pubr()

Figure1 <- ggarrange(F1, F2, F3, F4,
                     ncol = 4, nrow = 1, 
                     common.legend = TRUE)

annotate_figure(Figure1, 
                left = text_grob("Standardized Mean Differences\n for democratic breakdown/ recession", rot = 90))
ggsave("results/Figure_3/Figure_Covariate_Balance_Autocratization.pdf", height = 15, width = 35, units= c("cm"))


#### Improved Covariate Balance due to the Refinement of Matched Sets. 

par(mfrow=c(1,3)) 
balance_scatter(non_refined_set = PM.results.none$att,
                refined_list = list(PM.results_mahalanobis$att),
                data = vdem,
                covariates = c("e_gdppc",  "milper_extrapol", "eprratio_extrapol",
                               "e_pop",  "ma_count_attacks_dom", "democracy_stock", "e_area_log"))
balance_scatter(non_refined_set = PM.results.none$att,
                refined_list = list(PM.results_psmatch$att),
                data = vdem,
                covariates = c("e_gdppc",  "milper_extrapol", "eprratio_extrapol",
                               "e_pop",  "ma_count_attacks_dom", "democracy_stock", "e_area_log"))
balance_scatter(non_refined_set = PM.results.none$att,
                refined_list = list(PM.results_psweight$att),
                data = vdem,
                covariates = c("e_gdppc",  "milper_extrapol", "eprratio_extrapol",
                               "e_pop",  "ma_count_attacks_dom", "democracy_stock", "e_area_log"))

dev.copy(pdf,'results/Figure_3/Figure_Improved_Covariate_Balance.pdf', width = 10, height = 5)
dev.off()

##### Democratic Breakdown as Treatment and Terrorist Attacks ####

## domestic terrorist attacks ##

PM_M1_aut <- PanelMatch(lag = 2, time.id = "year", unit.id = "country_id", 
                        treatment = "democratic_breakdown_recession", refinement.method = "ps.weight", 
                        data = vdem, match.missing = FALSE, 
                        covs.formula = ~ I(lag(e_gdppc, lag_2)) + 
                          I(lag(milper_extrapol, lag_2)) + I(lag(eprratio_extrapol, lag_2)) + 
                          I(lag(e_pop, lag_2)) + I(lag(democracy_stock, lag_2)) + I(lag(e_area_log, lag_2)) + 
                          I(lag(ma_count_attacks_dom, lag_2)), 
                        size.match = 10, qoi = "att" ,outcome.var = "count_attacks_dom",
                        lead = lead_5, forbid.treatment.reversal = FALSE,
                        use.diagonal.variance.matrix = TRUE)

M1_aut<- PanelEstimate(sets = PM_M1_aut, data = vdem)

M1_aut_ci09 <- PanelEstimate(sets = PM_M1_aut, data = vdem, confidence.level = .9)


M1_placebo <- placebo_test(PM_M1_aut,
                           data = vdem,
                           lag.in = 2,
                           number.iterations = 1000,
                           confidence.level = 0.95, 
                           plot = F)

M1_placebo_t_2_est <- M1_placebo$estimates[1]
M1_placebo_t_2_high <- quantile(M1_placebo$bootstrapped.estimates[ ,1],0.975)
M1_placebo_t_2_low <- quantile(M1_placebo$bootstrapped.estimates[ ,1],0.025)

M1_placebo_t_2_high_90 <- quantile(M1_placebo$bootstrapped.estimates[ ,1],0.95)
M1_placebo_t_2_low_90 <- quantile(M1_placebo$bootstrapped.estimates[ ,1],0.05)


M1_placebo_data <- data.frame(year = c(-2), 
                              estimate =  M1_placebo_t_2_est,
                              conf_low_95 = M1_placebo_t_2_low,
                              conf_low_90 = M1_placebo_t_2_low_90,
                              conf_high_95 = M1_placebo_t_2_high, 
                              conf_high_90 = M1_placebo_t_2_high_90)


## terrorist fatalities domestic ##

PM_M2_aut <- PanelMatch(lag =2, time.id = "year", unit.id = "country_id", 
                        treatment = "democratic_breakdown_recession", refinement.method = "ps.weight", 
                        data = vdem, match.missing = FALSE, 
                        covs.formula = ~ I(lag(e_gdppc, lag_2)) + 
                          I(lag(milper_extrapol, lag_2)) + I(lag(eprratio_extrapol, lag_2)) + 
                          I(lag(e_pop, lag_2)) + I(lag(democracy_stock, lag_2)) + I(lag(e_area_log, lag_2)) + 
                          I(lag(ma_sum_nkill_dom, lag_2)), 
                        size.match = 10, qoi = "att" ,outcome.var = "sum_nkill_dom",
                        lead = lead_5, forbid.treatment.reversal = FALSE,
                        use.diagonal.variance.matrix = TRUE)
M2_aut <- PanelEstimate(sets = PM_M2_aut, data = vdem)
M2_aut_ci09 <- PanelEstimate(sets = PM_M2_aut, data = vdem, confidence.level = .9)

par(mfrow=c(1,2)) 
plot(M1_aut)
plot(M2_aut)

#### Save results as data.frames to build coherent predicted plots with 2WFE and PanelMatch Models ####

#95% CIs ##

M1_pred <- summary(M1_aut)

M1_pred <- M1_pred$summary

M1_pred <- cbind(M1_pred , time = c(0,1,2,3,4,5))

M2_pred <- summary(M2_aut)

M2_pred <- M2_pred$summary

M2_pred <- cbind(M2_pred , time = c(0,1,2,3,4,5))


saveRDS(M1_pred, "results/Figure_3/data_M1_pred_dem_breakdown.rds")
saveRDS(M2_pred, "results/Figure_3/data_M2_pred_dem_breakdown.rds")

# 90% CIs ##

M1_pred_ci09 <- summary(M1_aut_ci09)

M1_pred_ci09 <- M1_pred_ci09$summary

M1_pred_ci09 <- cbind(M1_pred_ci09 , time = c(0,1,2,3,4,5))

M2_pred_ci09 <- summary(M2_aut_ci09)

M2_pred_ci09 <- M2_pred_ci09$summary

M2_pred_ci09 <- cbind(M2_pred_ci09 , time = c(0,1,2,3,4,5))

saveRDS(M1_pred_ci09, "results/Figure_3/data_M1_pred_dem_breakdown_ci09.rds")
saveRDS(M2_pred_ci09, "results/Figure_3/data_M2_pred_dem_breakdown_ci09.rds")


#### Regressed Autocracy Effect ####

data_figure_1a <- readRDS("results/Figure_3/data_M1_pred_dem_breakdown.rds") 
data_figure_1b <- readRDS("results/Figure_3/data_M2_pred_dem_breakdown.rds") 

data_figure_1a <- as.data.frame(data_figure_1a)
data_figure_1b <- as.data.frame(data_figure_1b)

data_figure_1a_09 <- readRDS("results/Figure_3/data_M1_pred_dem_breakdown_ci09.rds") 
data_figure_1b_09 <- readRDS("results/Figure_3/data_M2_pred_dem_breakdown_ci09.rds") 

data_figure_1a_09 <- as.data.frame(data_figure_1a_09)
data_figure_1b_09 <- as.data.frame(data_figure_1b_09)

data_figure_1a <- cbind(data_figure_1a, data_figure_1a_09[,3:4])
data_figure_1b <- cbind(data_figure_1b, data_figure_1b_09[,3:4])


data_figure_1a <- data_figure_1a  %>%
  rename(SE = std.error, 
         conf_low_95 = `2.5%`, 
         conf_high_95 = `97.5%`,
         conf_low_90 = `5%`, 
         conf_high_90 = `95%`, 
         year = time) %>%
  dplyr::select(year, estimate, conf_low_95, conf_low_90, conf_high_95, conf_high_90)


data_figure_1a <- rbind(data_figure_1a, M1_placebo_data)

data_figure_1a[nrow(data_figure_1a) + 1,] <- c(-1, 0, NA,  NA, NA, NA)


#### Building Subplots for each Estimand ####

f1a <- ggplot(data = data_figure_1a, aes(x = year)) +
  geom_linerange(aes(ymin=conf_low_95, ymax = conf_high_95), size = 1.3, alpha = 0.7, color = "#a6a6a6") +
  geom_linerange(aes(ymin=conf_low_90, ymax = conf_high_90), size = 1.5, color = "#5a5a5a") +
  geom_point(aes(y=estimate), size = 4, color = "black") +
  geom_hline(aes(yintercept = 0), color = "red", linetype = "dashed", alpha = 0.5) +
  theme_pubr()+
  scale_x_continuous(breaks=seq(-2,10,1)) +
  labs(x = "Year around Democratic Breakdown/Recession", 
       y = "Estimated effect of Democratic Breakdown/Recession", 
       title = "Propensity Score Weighting: Terrorist Attacks")

f1a 

ggsave("outputs/Figure_3.png", units= c("cm"), width = 20, height = 12, dpi = 900)
ggsave("outputs/Figure_3.pdf", units= c("cm"), width = 20, height = 12, dpi = 900)

#------------------------------------------------------------------------------#
#---------------------------------- Figure 4 ----------------------------------#
#------------------------------------------------------------------------------#

vdem <- distinct(vdem, country_id, year, .keep_all= TRUE) # Distinct observations

vdem$year <- as.integer(vdem$year)


vdem <- as.data.frame(vdem)
class(vdem)

# Define Lags and leads for all

lag_2 <- c(1:2)
lead_5 <- c(0:5)

vdem$country_id <- as.integer(vdem$country_id)

summary(vdem$regressed_autocracy)
table(vdem$regressed_autocracy)

#### 1. Effect of Autocratization on Terrorism in all regimes #####

## Before Refinement ##
PM.results.none <- PanelMatch(lag = 2, time.id = "year", unit.id = "country_id", 
                              treatment = "regressed_autocracy", refinement.method = "none", 
                              data = vdem, match.missing = FALSE, 
                              covs.formula = ~ I(lag(ma_count_attacks_dom, lag_2)), 
                              size.match = 10, qoi = "att" ,outcome.var = "count_attacks_dom",
                              lead = lead_5, forbid.treatment.reversal = FALSE, 
                              use.diagonal.variance.matrix = TRUE)

## Mahanlanobis Distance Matching ##
PM.results_mahalanobis <- PanelMatch(lag = 2, time.id = "year", unit.id = "country_id", 
                                     treatment = "regressed_autocracy", refinement.method = "mahalanobis", 
                                     data = vdem, match.missing = FALSE, 
                                     covs.formula = ~ I(lag(e_gdppc, lag_2)) + 
                                       I(lag(milper_extrapol, lag_2)) + I(lag(eprratio_extrapol, lag_2)) + 
                                       I(lag(e_pop, lag_2)) + I(lag(democracy_stock, lag_2)) + I(lag(e_area_log, lag_2)) + 
                                       I(lag(ma_count_attacks_dom, lag_2)), 
                                     size.match = 10, qoi = "att" ,outcome.var = "count_attacks_dom",
                                     lead = lead_5, forbid.treatment.reversal = FALSE, 
                                     use.diagonal.variance.matrix = TRUE)

## Propensity Score Matching ##
PM.results_psmatch <- PanelMatch(lag = 2, time.id = "year", unit.id = "country_id", 
                                 treatment = "regressed_autocracy", refinement.method = "ps.match", 
                                 data = vdem, match.missing = FALSE, 
                                 covs.formula = ~ I(lag(e_gdppc, lag_2)) + 
                                   I(lag(milper_extrapol, lag_2)) + I(lag(eprratio_extrapol, lag_2)) + 
                                   I(lag(e_pop, lag_2)) + I(lag(democracy_stock, lag_2)) + I(lag(e_area_log, lag_2)) + 
                                   I(lag(ma_count_attacks_dom, lag_2)), 
                                 size.match = 10, qoi = "att" ,outcome.var = "count_attacks_dom",
                                 lead = lead_5, forbid.treatment.reversal = FALSE,
                                 use.diagonal.variance.matrix = TRUE)

## Propensity Score Weighting ##
PM.results_psweight <- PanelMatch(lag = 2, time.id = "year", unit.id = "country_id", 
                                  treatment = "regressed_autocracy", refinement.method = "ps.weight", 
                                  data = vdem, match.missing = FALSE, 
                                  covs.formula = ~ I(lag(e_gdppc, lag_2)) + 
                                    I(lag(milper_extrapol, lag_2)) + I(lag(eprratio_extrapol, lag_2)) + 
                                    I(lag(e_pop, lag_2)) + I(lag(democracy_stock, lag_2)) + I(lag(e_area_log, lag_2)) + 
                                    I(lag(ma_count_attacks_dom, lag_2)), 
                                  size.match = 10, qoi = "att" ,outcome.var = "count_attacks_dom",
                                  lead = lead_5, forbid.treatment.reversal = FALSE,
                                  use.diagonal.variance.matrix = TRUE)

class(PM.results_psweight)

msets <- PM.results_psweight$att

print(msets)

######################################################

## Display different matched sets ##
mset <- PM.results_psweight$att[10]

DisplayTreatment(unit.id = "country_id",
                 time.id = "year", legend.position = "none",
                 xlab = "year", ylab = "Country Code",
                 treatment = "regressed_autocracy", data = vdem,
                 matched.set = mset, # this way we highlight the particular set
                 show.set.only = TRUE)

summary(PM.results_mahalanobis$att)
summary(PM.results_psweight$att)


par(mfrow=c(1,1)) 
plot(PM.results_psweight$att, main ="Distribution of Matched Set Sizes")

dev.copy(pdf,'results/Figure_4/Figure_Frequency_Distribution.pdf', width = 10, height = 5)
dev.off()


#### Check Improved Covariate Balance due to Matching over the Pre-Treatment Time ####

## Get covariate balance ##

df_none <- get_covariate_balance(PM.results.none$att,
                                 data = vdem,
                                 covariates = c("e_gdppc",  "milper_extrapol", "eprratio_extrapol",
                                                "e_pop", "ma_count_attacks_dom", "democracy_stock", "e_area_log"),
                                 plot = FALSE)

df_mahalanobis <- get_covariate_balance(PM.results_mahalanobis$att,
                                        data = vdem,
                                        covariates = c("e_gdppc",  "milper_extrapol", "eprratio_extrapol",
                                                       "e_pop", "ma_count_attacks_dom", "democracy_stock", "e_area_log"),
                                        plot = FALSE)


df_psmatch <- get_covariate_balance(PM.results_psmatch$att,
                                    data = vdem,
                                    covariates = c("e_gdppc",  "milper_extrapol", "eprratio_extrapol",
                                                   "e_pop", "ma_count_attacks_dom", "democracy_stock", "e_area_log"),
                                    plot = FALSE)


df_psweight <- get_covariate_balance(PM.results_psweight$att,
                                     data = vdem,
                                     covariates = c("e_gdppc",  "milper_extrapol", "eprratio_extrapol",
                                                    "e_pop", "ma_count_attacks_dom", "democracy_stock", "e_area_log"),
                                     plot = FALSE)


#Tidying outputs #
df_none <- as.data.frame(df_none)

df_none <- df_none %>%
  rownames_to_column("year") %>%
  pivot_longer(-year, names_to = "covariate", values_to = "sd_covariate") %>%
  mutate(year = case_when(year =="t_4" ~ -4, 
                          year =="t_3" ~ -3,
                          year =="t_2" ~ -2,
                          year =="t_1" ~ -1,
                          year == "t_0" ~ 0))

F1 <- ggplot(df_none, aes(x = year, y = sd_covariate, color = covariate)) +
  geom_line() +
  geom_hline(yintercept=0, color = "black") +
  labs(x = "", y = "", 
       title = "Before Refinement") +
  ylim(-0.8, 0.8) +
  theme(legend.position = "none", 
        plot.title = element_text(hjust = 0.5)) +
  theme_pubr()

df_mahalanobis <- as.data.frame(df_mahalanobis)

df_mahalanobis <- df_mahalanobis %>%
  rownames_to_column("year") %>%
  pivot_longer(-year, names_to = "covariate", values_to = "sd_covariate") %>%
  mutate(year = case_when(year =="t_4" ~ -4, 
                          year =="t_3" ~ -3,
                          year =="t_2" ~ -2,
                          year =="t_1" ~ -1,
                          year == "t_0" ~ 0))

F2 <- ggplot(df_mahalanobis, aes(x = year, y = sd_covariate, color = covariate)) +
  geom_line() +
  geom_hline(yintercept=0, color = "black") +
  labs(x = "", y = "", 
       title = "Mahanlanobis Distance\n Matching") +
  ylim(-0.8, 0.8) +
  theme(legend.position = "none", 
        plot.title = element_text(hjust = 0.5)) +
  theme_pubr()

df_psmatch <- as.data.frame(df_psmatch)

df_psmatch <- df_psmatch %>%
  rownames_to_column("year") %>%
  pivot_longer(-year, names_to = "covariate", values_to = "sd_covariate") %>%
  mutate(year = case_when(year =="t_4" ~ -4, 
                          year =="t_3" ~ -3,
                          year =="t_2" ~ -2,
                          year =="t_1" ~ -1,
                          year == "t_0" ~ 0))

F3 <- ggplot(df_psmatch, aes(x = year, y = sd_covariate, color = covariate)) +
  geom_line() +
  geom_hline(yintercept=0, color = "black") +
  labs(x = "", y = "", 
       title = "Propensity Score\n Matching") +
  ylim(-0.8, 0.8) +
  theme(legend.position = "none", 
        plot.title = element_text(hjust = 0.5)) +
  theme_pubr()


df_psweight <- as.data.frame(df_psweight)

df_psweight <- df_psweight %>%
  rownames_to_column("year") %>%
  pivot_longer(-year, names_to = "covariate", values_to = "sd_covariate") %>%
  mutate(year = case_when(year =="t_4" ~ -4, 
                          year =="t_3" ~ -3,
                          year =="t_2" ~ -2,
                          year =="t_1" ~ -1,
                          year == "t_0" ~ 0))

F4 <- ggplot(df_psweight, aes(x = year, y = sd_covariate, color = covariate)) +
  geom_line() +
  geom_hline(yintercept=0, color = "black") +
  labs(x = "", y = "", 
       title = "Propensity Score\n Weighting") +
  ylim(-0.8, 0.8) +
  theme(legend.position = "none", 
        plot.title = element_text(hjust = 0.5)) +
  theme_pubr()

Figure1 <- ggarrange(F1, F2, F3, F4,
                     ncol = 4, nrow = 1, 
                     common.legend = TRUE)

annotate_figure(Figure1, 
                left = text_grob("Standardized Mean Differences\n for Autocratization", rot = 90))
ggsave("results/Figure_4/Figure_Covariate_Balance_Autocratization.pdf", height = 15, width = 35, units= c("cm"))


#### Improved Covariate Balance due to the Refinement of Matched Sets. 

par(mfrow=c(1,3)) 
balance_scatter(non_refined_set = PM.results.none$att,
                refined_list = list(PM.results_mahalanobis$att),
                data = vdem,
                covariates = c("e_gdppc",  "milper_extrapol", "eprratio_extrapol",
                               "e_pop",  "ma_count_attacks_dom", "democracy_stock", "e_area_log"))
balance_scatter(non_refined_set = PM.results.none$att,
                refined_list = list(PM.results_psmatch$att),
                data = vdem,
                covariates = c("e_gdppc",  "milper_extrapol", "eprratio_extrapol",
                               "e_pop",  "ma_count_attacks_dom", "democracy_stock", "e_area_log"))
balance_scatter(non_refined_set = PM.results.none$att,
                refined_list = list(PM.results_psweight$att),
                data = vdem,
                covariates = c("e_gdppc",  "milper_extrapol", "eprratio_extrapol",
                               "e_pop",  "ma_count_attacks_dom", "democracy_stock", "e_area_log"))

dev.copy(pdf,'results/Figure_4/Figure_Improved_Covariate_Balance.pdf', width = 10, height = 5)
dev.off()

##### Democratic Breakdown as Treatment and Terrorist Attacks ####

## domestic terrorist attacks ##

PM_M1_aut <- PanelMatch(lag = 2, time.id = "year", unit.id = "country_id", 
                        treatment = "regressed_autocracy", refinement.method = "ps.weight", 
                        data = vdem, match.missing = FALSE, 
                        covs.formula = ~ I(lag(e_gdppc, lag_2)) + 
                          I(lag(milper_extrapol, lag_2)) + I(lag(eprratio_extrapol, lag_2)) + 
                          I(lag(e_pop, lag_2)) + I(lag(democracy_stock, lag_2)) + I(lag(e_area_log, lag_2)) + 
                          I(lag(ma_count_attacks_dom, lag_2)), 
                        size.match = 10, qoi = "att" ,outcome.var = "count_attacks_dom",
                        lead = lead_5, forbid.treatment.reversal = FALSE,
                        use.diagonal.variance.matrix = TRUE)

M1_aut<- PanelEstimate(sets = PM_M1_aut, data = vdem)

M1_aut_ci09 <- PanelEstimate(sets = PM_M1_aut, data = vdem, confidence.level = .9)


M1_placebo <- placebo_test(PM_M1_aut,
                           data = vdem,
                           lag.in = 2,
                           number.iterations = 1000,
                           confidence.level = 0.95, 
                           plot = F)

M1_placebo_t_2_est <- M1_placebo$estimates[1]
M1_placebo_t_2_high <- quantile(M1_placebo$bootstrapped.estimates[ ,1],0.975)
M1_placebo_t_2_low <- quantile(M1_placebo$bootstrapped.estimates[ ,1],0.025)

M1_placebo_t_2_high_90 <- quantile(M1_placebo$bootstrapped.estimates[ ,1],0.95)
M1_placebo_t_2_low_90 <- quantile(M1_placebo$bootstrapped.estimates[ ,1],0.05)


M1_placebo_data <- data.frame(year = c(-2), 
                              estimate =  M1_placebo_t_2_est,
                              conf_low_95 = M1_placebo_t_2_low,
                              conf_low_90 = M1_placebo_t_2_low_90,
                              conf_high_95 = M1_placebo_t_2_high, 
                              conf_high_90 = M1_placebo_t_2_high_90)


## terrorist fatalities domestic ##

PM_M2_aut <- PanelMatch(lag =2, time.id = "year", unit.id = "country_id", 
                        treatment = "regressed_autocracy", refinement.method = "ps.weight", 
                        data = vdem, match.missing = FALSE, 
                        covs.formula = ~ I(lag(e_gdppc, lag_2)) + 
                          I(lag(milper_extrapol, lag_2)) + I(lag(eprratio_extrapol, lag_2)) + 
                          I(lag(e_pop, lag_2)) + I(lag(democracy_stock, lag_2)) + I(lag(e_area_log, lag_2)) + 
                          I(lag(ma_sum_nkill_dom, lag_2)), 
                        size.match = 10, qoi = "att" ,outcome.var = "sum_nkill_dom",
                        lead = lead_5, forbid.treatment.reversal = FALSE,
                        use.diagonal.variance.matrix = TRUE)
M2_aut <- PanelEstimate(sets = PM_M2_aut, data = vdem)
M2_aut_ci09 <- PanelEstimate(sets = PM_M2_aut, data = vdem, confidence.level = .9)

par(mfrow=c(1,2)) 
plot(M1_aut)
plot(M2_aut)

#### Save results as data.frames to build coherent predicted plots with 2WFE and PanelMatch Models ####

#95% CIs ##

M1_pred <- summary(M1_aut)

M1_pred <- M1_pred$summary

M1_pred <- cbind(M1_pred , time = c(0,1,2,3,4,5))

M2_pred <- summary(M2_aut)

M2_pred <- M2_pred$summary

M2_pred <- cbind(M2_pred , time = c(0,1,2,3,4,5))


saveRDS(M1_pred, "results/Figure_4/data_M1_pred_dem_breakdown.rds")
saveRDS(M2_pred, "results/Figure_4/data_M2_pred_dem_breakdown.rds")

# 90% CIs ##

M1_pred_ci09 <- summary(M1_aut_ci09)

M1_pred_ci09 <- M1_pred_ci09$summary

M1_pred_ci09 <- cbind(M1_pred_ci09 , time = c(0,1,2,3,4,5))

M2_pred_ci09 <- summary(M2_aut_ci09)

M2_pred_ci09 <- M2_pred_ci09$summary

M2_pred_ci09 <- cbind(M2_pred_ci09 , time = c(0,1,2,3,4,5))

saveRDS(M1_pred_ci09, "results/Figure_4/data_M1_pred_dem_breakdown_ci09.rds")
saveRDS(M2_pred_ci09, "results/Figure_4/data_M2_pred_dem_breakdown_ci09.rds")


#### Regressed Autocracy Effect ####

data_figure_1a <- readRDS("results/Figure_4/data_M1_pred_dem_breakdown.rds") 
data_figure_1b <- readRDS("results/Figure_4/data_M2_pred_dem_breakdown.rds") 

data_figure_1a <- as.data.frame(data_figure_1a)
data_figure_1b <- as.data.frame(data_figure_1b)

data_figure_1a_09 <- readRDS("results/Figure_4/data_M1_pred_dem_breakdown_ci09.rds") 
data_figure_1b_09 <- readRDS("results/Figure_4/data_M2_pred_dem_breakdown_ci09.rds") 

data_figure_1a_09 <- as.data.frame(data_figure_1a_09)
data_figure_1b_09 <- as.data.frame(data_figure_1b_09)

data_figure_1a <- cbind(data_figure_1a, data_figure_1a_09[,3:4])
data_figure_1b <- cbind(data_figure_1b, data_figure_1b_09[,3:4])

data_figure_1a <- data_figure_1a  %>%
  rename(SE = std.error, 
         conf_low_95 = `2.5%`, 
         conf_high_95 = `97.5%`,
         conf_low_90 = `5%`, 
         conf_high_90 = `95%`, 
         year = time) %>%
  dplyr::select(year, estimate, conf_low_95, conf_low_90, conf_high_95, conf_high_90)


data_figure_1a <- rbind(data_figure_1a, M1_placebo_data)

data_figure_1a[nrow(data_figure_1a) + 1,] <- c(-1, 0, NA,  NA, NA, NA)


#### Building Subplots for each Estimand ####

f1a <- ggplot(data = data_figure_1a, aes(x = year)) +
  geom_linerange(aes(ymin=conf_low_95, ymax = conf_high_95), size = 1.3, alpha = 0.7, color = "#a6a6a6") +
  geom_linerange(aes(ymin=conf_low_90, ymax = conf_high_90), size = 1.5, color = "#5a5a5a") +
  geom_point(aes(y=estimate), size = 4, color = "black") +
  geom_hline(aes(yintercept = 0), color = "red", linetype = "dashed", alpha = 0.5) +
  theme_pubr()+
  scale_x_continuous(breaks=seq(-2,5,1)) +
  labs(x = "Year around Autocratic Consolidation", 
       y = "Estimated effect of Autocratic Consolidation", 
       title = "Propensity Score Weighting: Terrorist Attacks")
f1a

ggsave("outputs/Figure_4.png", units= c("cm"), width = 20, height = 12, dpi = 900)
ggsave("outputs/Figure_4.pdf", units= c("cm"), width = 20, height = 12, dpi = 900)


#------------------------------------------------------------------------------#
#---------------------------------- Figure 5 ----------------------------------#
#------------------------------------------------------------------------------#

vdem <- distinct(vdem, country_id, year, .keep_all= TRUE) # Distinct observations

vdem$year <- as.integer(vdem$year)


vdem <- as.data.frame(vdem)
class(vdem)



test <- vdem %>%
  filter(year >=1970) %>%
  filter(year <=2020) %>%
  filter(aut_ep_id !=0) %>%
  group_by(aut_ep_id) %>%
  summarize(v2x_regime = first(v2x_regime)) %>%
  drop_na(aut_ep_id)

table(test$v2x_regime)

#### descriptives ####

vdem_country_examples <- vdem %>%
  filter(year >=1970) %>%
  dplyr::select(country_name, aut_ep, aut_ep_id, v2x_regime, count_attacks_dom, everything()) %>%
  mutate(v2x_regime = as.factor(v2x_regime)) %>%
  arrange(v2x_regime, aut_ep_id, desc(count_attacks_dom))


# Define Lags and leads for all

lag_2 <- c(1:2)
lead_5 <- c(0:5)

vdem <- vdem %>%
  dplyr::select(-c(country_name, aut_ep_id))

vdem$country_id <- as.integer(vdem$country_id)

#### 1. Effect of democratic_breakdown_delta on Terrorism in all regimes #####

## Before Refinement ##
PM.results.none <- PanelMatch(lag = 2, time.id = "year", unit.id = "country_id", 
                              treatment = "autocratization_treatment", refinement.method = "none", 
                              data = vdem, match.missing = FALSE, 
                              covs.formula = ~ I(lag(ma_count_attacks_dom, lag_2)), 
                              size.match = 10, qoi = "att" ,outcome.var = "count_attacks_dom",
                              lead = lead_5, forbid.treatment.reversal = FALSE, 
                              use.diagonal.variance.matrix = TRUE)


## Mahanlanobis Distance Matching ##
PM.results_mahalanobis <- PanelMatch(lag = 2, time.id = "year", unit.id = "country_id", 
                                     treatment = "autocratization_treatment", refinement.method = "mahalanobis", 
                                     data = vdem, match.missing = FALSE, 
                                     covs.formula = ~ I(lag(e_gdppc, lag_2)) + 
                                       I(lag(milper_extrapol, lag_2)) + I(lag(eprratio_extrapol, lag_2)) + 
                                       I(lag(e_pop, lag_2)) + I(lag(democracy_stock, lag_2)) + I(lag(e_area_log, lag_2)) + 
                                       I(lag(ma_count_attacks_dom, lag_2)), 
                                     size.match = 10, qoi = "att" ,outcome.var = "count_attacks_dom",
                                     lead = lead_5, forbid.treatment.reversal = FALSE, 
                                     use.diagonal.variance.matrix = TRUE)

## Propensity Score Matching ##
PM.results_psmatch <- PanelMatch(lag = 2, time.id = "year", unit.id = "country_id", 
                                 treatment = "autocratization_treatment", refinement.method = "ps.match", 
                                 data = vdem, match.missing = FALSE, 
                                 covs.formula = ~ I(lag(e_gdppc, lag_2)) + 
                                   I(lag(milper_extrapol, lag_2)) + I(lag(eprratio_extrapol, lag_2)) + 
                                   I(lag(e_pop, lag_2)) + I(lag(democracy_stock, lag_2)) + I(lag(e_area_log, lag_2)) + 
                                   I(lag(ma_count_attacks_dom, lag_2)), 
                                 size.match = 10, qoi = "att" ,outcome.var = "count_attacks_dom",
                                 lead = lead_5, forbid.treatment.reversal = FALSE,
                                 use.diagonal.variance.matrix = TRUE)

## Propensity Score Weighting ##
PM.results_psweight <- PanelMatch(lag = 2, time.id = "year", unit.id = "country_id", 
                                  treatment = "autocratization_treatment", refinement.method = "ps.weight", 
                                  data = vdem, match.missing = FALSE, 
                                  covs.formula = ~ I(lag(e_gdppc, lag_2)) + 
                                    I(lag(milper_extrapol, lag_2)) + I(lag(eprratio_extrapol, lag_2)) + 
                                    I(lag(e_pop, lag_2)) + I(lag(democracy_stock, lag_2)) + I(lag(e_area_log, lag_2)) + 
                                    I(lag(ma_count_attacks_dom, lag_2)), 
                                  size.match = 10, qoi = "att" ,outcome.var = "count_attacks_dom",
                                  lead = lead_5, forbid.treatment.reversal = FALSE,
                                  use.diagonal.variance.matrix = TRUE)

class(PM.results_psweight)

msets <- PM.results_psweight$att

print(msets)

######################################################

## Display different matched sets ##
mset <- PM.results_psweight$att[10]

DisplayTreatment(unit.id = "country_id",
                 time.id = "year", legend.position = "none",
                 xlab = "year", ylab = "Country Code",
                 treatment = "autocratization_treatment", data = vdem,
                 matched.set = mset, # this way we highlight the particular set
                 show.set.only = TRUE)

summary(PM.results_mahalanobis$att)
summary(PM.results_psweight$att)

par(mfrow=c(1,1)) 
plot(PM.results_psweight$att, main ="Distribution of Matched Set Sizes")

dev.copy(pdf,'results/Figure_5/Figure_Frequency_Distribution.pdf', width = 10, height = 5)
dev.off()


#### Check Improved Covariate Balance due to Matching over the Pre-Treatment Time ####

## Get covariate balance ##

df_none <- get_covariate_balance(PM.results.none$att,
                                 data = vdem,
                                 covariates = c("e_gdppc",  "milper_extrapol", "eprratio_extrapol",
                                                "e_pop", "ma_count_attacks_dom", "democracy_stock", "e_area_log"),
                                 plot = FALSE)

df_mahalanobis <- get_covariate_balance(PM.results_mahalanobis$att,
                                        data = vdem,
                                        covariates = c("e_gdppc",  "milper_extrapol", "eprratio_extrapol",
                                                       "e_pop", "ma_count_attacks_dom", "democracy_stock", "e_area_log"),
                                        plot = FALSE)


df_psmatch <- get_covariate_balance(PM.results_psmatch$att,
                                    data = vdem,
                                    covariates = c("e_gdppc",  "milper_extrapol", "eprratio_extrapol",
                                                   "e_pop", "ma_count_attacks_dom", "democracy_stock", "e_area_log"),
                                    plot = FALSE)


df_psweight <- get_covariate_balance(PM.results_psweight$att,
                                     data = vdem,
                                     covariates = c("e_gdppc",  "milper_extrapol", "eprratio_extrapol",
                                                    "e_pop", "ma_count_attacks_dom", "democracy_stock", "e_area_log"),
                                     plot = FALSE)


#Tidying outputs #
df_none <- as.data.frame(df_none)

df_none <- df_none %>%
  rownames_to_column("year") %>%
  pivot_longer(-year, names_to = "covariate", values_to = "sd_covariate") %>%
  mutate(year = case_when(year =="t_4" ~ -4, 
                          year =="t_3" ~ -3,
                          year =="t_2" ~ -2,
                          year =="t_1" ~ -1,
                          year == "t_0" ~ 0))

F1 <- ggplot(df_none, aes(x = year, y = sd_covariate, color = covariate)) +
  geom_line() +
  geom_hline(yintercept=0, color = "black") +
  labs(x = "", y = "", 
       title = "Before Refinement") +
  ylim(-0.8, 0.8) +
  theme(legend.position = "none", 
        plot.title = element_text(hjust = 0.5)) +
  theme_pubr()

df_mahalanobis <- as.data.frame(df_mahalanobis)

df_mahalanobis <- df_mahalanobis %>%
  rownames_to_column("year") %>%
  pivot_longer(-year, names_to = "covariate", values_to = "sd_covariate") %>%
  mutate(year = case_when(year =="t_4" ~ -4, 
                          year =="t_3" ~ -3,
                          year =="t_2" ~ -2,
                          year =="t_1" ~ -1,
                          year == "t_0" ~ 0))

F2 <- ggplot(df_mahalanobis, aes(x = year, y = sd_covariate, color = covariate)) +
  geom_line() +
  geom_hline(yintercept=0, color = "black") +
  labs(x = "", y = "", 
       title = "Mahanlanobis Distance\n Matching") +
  ylim(-0.8, 0.8) +
  theme(legend.position = "none", 
        plot.title = element_text(hjust = 0.5)) +
  theme_pubr()

df_psmatch <- as.data.frame(df_psmatch)

df_psmatch <- df_psmatch %>%
  rownames_to_column("year") %>%
  pivot_longer(-year, names_to = "covariate", values_to = "sd_covariate") %>%
  mutate(year = case_when(year =="t_4" ~ -4, 
                          year =="t_3" ~ -3,
                          year =="t_2" ~ -2,
                          year =="t_1" ~ -1,
                          year == "t_0" ~ 0))

F3 <- ggplot(df_psmatch, aes(x = year, y = sd_covariate, color = covariate)) +
  geom_line() +
  geom_hline(yintercept=0, color = "black") +
  labs(x = "", y = "", 
       title = "Propensity Score\n Matching") +
  ylim(-0.8, 0.8) +
  theme(legend.position = "none", 
        plot.title = element_text(hjust = 0.5)) +
  theme_pubr()


df_psweight <- as.data.frame(df_psweight)

df_psweight <- df_psweight %>%
  rownames_to_column("year") %>%
  pivot_longer(-year, names_to = "covariate", values_to = "sd_covariate") %>%
  mutate(year = case_when(year =="t_4" ~ -4, 
                          year =="t_3" ~ -3,
                          year =="t_2" ~ -2,
                          year =="t_1" ~ -1,
                          year == "t_0" ~ 0))

F4 <- ggplot(df_psweight, aes(x = year, y = sd_covariate, color = covariate)) +
  geom_line() +
  geom_hline(yintercept=0, color = "black") +
  labs(x = "", y = "", 
       title = "Propensity Score\n Weighting") +
  ylim(-0.8, 0.8) +
  theme(legend.position = "none", 
        plot.title = element_text(hjust = 0.5)) +
  theme_pubr()

Figure1 <- ggarrange(F1, F2, F3, F4,
                     ncol = 4, nrow = 1, 
                     common.legend = TRUE)

annotate_figure(Figure1, 
                left = text_grob("Standardized Mean Differences\n for autocratization", rot = 90))
ggsave("results/Figure_5/Figure_Covariate_Balance_Autocratization.pdf", height = 15, width = 35, units= c("cm"))


#### Improved Covariate Balance due to the Refinement of Matched Sets. 

par(mfrow=c(1,3)) 
balance_scatter(non_refined_set = PM.results.none$att,
                refined_list = list(PM.results_mahalanobis$att),
                data = vdem,
                covariates = c("e_gdppc",  "milper_extrapol", "eprratio_extrapol",
                               "e_pop",  "ma_count_attacks_dom",  "democracy_stock", "e_area_log"))
balance_scatter(non_refined_set = PM.results.none$att,
                refined_list = list(PM.results_psmatch$att),
                data = vdem,
                covariates = c("e_gdppc",  "milper_extrapol", "eprratio_extrapol",
                               "e_pop",  "ma_count_attacks_dom",  "democracy_stock", "e_area_log"))
balance_scatter(non_refined_set = PM.results.none$att,
                refined_list = list(PM.results_psweight$att),
                data = vdem,
                covariates = c("e_gdppc",  "milper_extrapol", "eprratio_extrapol",
                               "e_pop",  "ma_count_attacks_dom",  "democracy_stock", "e_area_log"))

dev.copy(pdf,'results/Figure_5/Figure_Improved_Covariate_Balance.pdf', width = 10, height = 5)
dev.off()

##### Democratic Breakdown as Treatment and Terrorist Attacks ####

## domestic terrorist attacks ##

PM_M1_aut <- PanelMatch::PanelMatch(lag = 2, time.id = "year", unit.id = "country_id", 
                        treatment = "autocratization_treatment", refinement.method = "ps.weight", 
                        data = vdem, match.missing = FALSE, 
                        covs.formula = ~ I(lag(e_gdppc, lag_2)) + 
                          I(lag(milper_extrapol, lag_2)) + I(lag(eprratio_extrapol, lag_2)) + 
                          I(lag(e_pop, lag_2)) +   I(lag(democracy_stock, lag_2)) + 
                          I(lag(e_area_log, lag_2)) + 
                          I(lag(ma_count_attacks_dom, lag_2)), 
                        size.match = 10, qoi = "att" ,outcome.var = "count_attacks_dom",
                        lead = lead_5, forbid.treatment.reversal = FALSE,
                        use.diagonal.variance.matrix = TRUE, placebo.test = T)

M1_aut<- PanelEstimate(sets = PM_M1_aut, data = vdem)

M1_aut_ci09 <- PanelEstimate(sets = PM_M1_aut, data = vdem, confidence.level = .9)

M1_placebo <- placebo_test(PM_M1_aut,
             data = vdem,
             lag.in = 2,
             number.iterations = 1000,
             confidence.level = 0.95, 
             plot = F)

M1_placebo_t_2_est <- M1_placebo$estimates[1]
M1_placebo_t_2_high <- quantile(M1_placebo$bootstrapped.estimates[ ,1],0.975)
M1_placebo_t_2_low <- quantile(M1_placebo$bootstrapped.estimates[ ,1],0.025)

M1_placebo_t_2_high_90 <- quantile(M1_placebo$bootstrapped.estimates[ ,1],0.95)
M1_placebo_t_2_low_90 <- quantile(M1_placebo$bootstrapped.estimates[ ,1],0.05)


M1_placebo_data <- data.frame(year = c(-2), 
                        estimate =  M1_placebo_t_2_est,
                        conf_low_95 = M1_placebo_t_2_low,
                        conf_low_90 = M1_placebo_t_2_low_90,
                        conf_high_95 = M1_placebo_t_2_high, 
                        conf_high_90 = M1_placebo_t_2_high_90)

## terrorist fatalities domestic ##

PM_M2_aut <- PanelMatch(lag =2, time.id = "year", unit.id = "country_id", 
                        treatment = "autocratization_treatment", refinement.method = "ps.weight", 
                        data = vdem, match.missing = FALSE, 
                        covs.formula = ~ I(lag(e_gdppc, lag_2)) + 
                          I(lag(milper_extrapol, lag_2)) + I(lag(eprratio_extrapol, lag_2)) + 
                          I(lag(e_pop, lag_2)) +   I(lag(democracy_stock, lag_2)) + I(lag(e_area_log, lag_2)) + 
                          I(lag(ma_sum_nkill_dom, lag_2)), 
                        size.match = 10, qoi = "att" ,outcome.var = "sum_nkill_dom",
                        lead = lead_5, forbid.treatment.reversal = FALSE,
                        use.diagonal.variance.matrix = TRUE)
M2_aut <- PanelEstimate(sets = PM_M2_aut, data = vdem)
M2_aut_ci09 <- PanelEstimate(sets = PM_M2_aut, data = vdem, confidence.level = .9)

par(mfrow=c(1,2)) 
plot(M1_aut)
plot(M2_aut)

#### Save results as data.frames to build coherent predicted plots with 2WFE and PanelMatch Models ####

#95% CIs ##

M1_pred <- summary(M1_aut)

M1_pred <- M1_pred$summary

M1_pred <- cbind(M1_pred , time = c(0,1,2,3,4,5))

M2_pred <- summary(M2_aut)

M2_pred <- M2_pred$summary

M2_pred <- cbind(M2_pred , time = c(0,1,2,3,4,5))


saveRDS(M1_pred, "results/Figure_5/data_M1_pred_dem_breakdown.rds")
saveRDS(M2_pred, "results/Figure_5/data_M2_pred_dem_breakdown.rds")

# 90% CIs ##

M1_pred_ci09 <- summary(M1_aut_ci09)

M1_pred_ci09 <- M1_pred_ci09$summary

M1_pred_ci09 <- cbind(M1_pred_ci09 , time = c(0,1,2,3,4,5))

M2_pred_ci09 <- summary(M2_aut_ci09)

M2_pred_ci09 <- M2_pred_ci09$summary

M2_pred_ci09 <- cbind(M2_pred_ci09 , time = c(0,1,2,3,4,5))

saveRDS(M1_pred_ci09, "results/Figure_5/data_M1_pred_dem_breakdown_ci09.rds")
saveRDS(M2_pred_ci09, "results/Figure_5/data_M2_pred_dem_breakdown_ci09.rds")


#### Autocratization Effect ####

data_figure_1a <- readRDS("results/Figure_5/data_M1_pred_dem_breakdown.rds") 
data_figure_1b <- readRDS("results/Figure_5/data_M2_pred_dem_breakdown.rds") 

data_figure_1a <- as.data.frame(data_figure_1a)
data_figure_1b <- as.data.frame(data_figure_1b)

data_figure_1a_09 <- readRDS("results/Figure_5/data_M1_pred_dem_breakdown_ci09.rds") 
data_figure_1b_09 <- readRDS("results/Figure_5/data_M2_pred_dem_breakdown_ci09.rds") 

data_figure_1a_09 <- as.data.frame(data_figure_1a_09)
data_figure_1b_09 <- as.data.frame(data_figure_1b_09)

data_figure_1a <- cbind(data_figure_1a, data_figure_1a_09[,3:4])
data_figure_1b <- cbind(data_figure_1b, data_figure_1b_09[,3:4])

data_figure_1a <- data_figure_1a  %>%
  rename(SE = std.error, 
         conf_low_95 = `2.5%`, 
         conf_high_95 = `97.5%`,
         conf_low_90 = `5%`, 
         conf_high_90 = `95%`, 
         year = time) %>%
  dplyr::select(year, estimate, conf_low_95, conf_low_90, conf_high_95, conf_high_90)
  
  
data_figure_1a <- rbind(data_figure_1a, M1_placebo_data)

data_figure_1a[nrow(data_figure_1a) + 1,] <- c(-1, 0, NA,  NA, NA, NA)



#### Building Subplots for each Estimand ####

f1a <- ggplot(data = data_figure_1a, aes(x = year)) +
  geom_linerange(aes(ymin=conf_low_95, ymax = conf_high_95), size = 1.3, alpha = 0.7, color = "#a6a6a6") +
  geom_linerange(aes(ymin= conf_low_90, ymax = conf_high_90), size = 1.5, color = "#5a5a5a") +
  geom_point(aes(y=estimate), size = 4, color = "black") +
  geom_hline(aes(yintercept = 0), color = "red", linetype = "dashed", alpha = 0.5) +
  theme_pubr()+
  scale_x_continuous(breaks=seq(-2,10,1)) +
  labs(x = "Year around Autocratization", 
       y = "Estimated effect of Autocratization", 
       title = "Propensity Score Weighting: Terrorist Attacks")

f1a
ggsave("outputs/Figure_5.png", units= c("cm"), width = 20, height = 12, dpi = 900)
ggsave("outputs/Figure_5.pdf", units= c("cm"), width = 20, height = 12, dpi = 900)



